


% HHVR_EM1b_3_types.m


  function  [loglik  ] = HHVR_EM1b_3_types(bbeta)  
  global DATA  X choices y_short zz individuals  

pi_u = 1/3; pi_e = 1/3;pi_f = 1/3; 
s=size(DATA,1)/choices;     % individuals x choices
nume1=zeros(s,1);nume2=zeros(s,1);nume3=zeros(s,1) ;nume4=zeros(s,1) ;nume5=zeros(s,1) ; 
numu1=zeros(s,1);numu2=zeros(s,1);numu3=zeros(s,1) ;numu4=zeros(s,1) ;numu5=zeros(s,1) ;
numf1=zeros(s,1);numf2=zeros(s,1);numf3=zeros(s,1) ;numf4=zeros(s,1) ;numf5=zeros(s,1) ;
prob_e1=zeros(s,1);prob_e2=zeros(s,1);prob_e3=zeros(s,1); prob_e4=zeros(s,1);prob_e5=zeros(s,1); 
prob_u1=zeros(s,1);prob_u2=zeros(s,1);prob_u3=zeros(s,1); prob_u4=zeros(s,1);prob_u5=zeros(s,1); 
prob_f1=zeros(s,1);prob_f2=zeros(s,1);prob_f3=zeros(s,1); prob_f4=zeros(s,1);prob_f5=zeros(s,1); 

denu= zeros(s,1); dene= zeros(s,1);   denf= zeros(s,1); 
denu_x1= zeros(s,1); dene_x1= zeros(s,1); denf_x1= zeros(s,1);

ze=zz(:,1); zu=zz(:,2);zf=zz(:,3); 
 
for i=1:(size(DATA,1))/choices                    % ie from 1 to number of individuals
       
     nume1 (i) =  exp(X(i*choices-4,1)*bbeta(5)    ) ;
     nume2 (i) =  exp(X(i*choices-3,1)*bbeta(5)  + 1*bbeta(1)  ) ;
     nume3 (i) =  exp(X(i*choices-2,1)*bbeta(5)  + 1*bbeta(2)  ) ;
     nume4 (i) =  exp(X(i*choices-1,1)*bbeta(5)  + 1*bbeta(3)  ) ;
     nume5 (i) =  exp(X(i*choices  ,1)*bbeta(5)  + 1*bbeta(4)  ) ;
     dene(i)   =  nume1 (i)  +  nume2 (i) + nume3 (i)+ nume4 (i)+ nume5 (i) ;
     dene_x1(i)   =  nume1(i)*X(i*choices-4,1)  +  nume2(i)*X(i*choices-3,1) + nume3(i)*X(i*choices-2,1)+ ... 
                     nume4(i)*X(i*choices-1,1) +   nume5(i)*X(i*choices  ,1 ) ;
     prob_e1(i)     =  nume1 (i) /  dene(i)  ;     prob_e2(i)     =  nume2 (i) /  dene(i)  ;
     prob_e3(i)     =  nume3 (i) /  dene(i)  ;     prob_e4(i)     =  nume4 (i) /  dene(i)  ;
     prob_e5(i)     =  nume5 (i) /  dene(i)  ;  
    
     
     numu1 (i) =  exp(X(i*choices-4,2)*bbeta(6)    ) ;
     numu2 (i) =  exp(X(i*choices-3,2)*bbeta(6)  + 1*bbeta(1)  ) ;
     numu3 (i) =  exp(X(i*choices-2,2)*bbeta(6)  + 1*bbeta(2)  ) ;
     numu4 (i) =  exp(X(i*choices-1,2)*bbeta(6)  + 1*bbeta(3)  ) ;
     numu5 (i) =  exp(X(i*choices  ,2)*bbeta(6)  + 1*bbeta(4)  ) ;
     denu(i)   =  numu1 (i)  +  numu2 (i) + numu3 (i)+ numu4 (i)+ numu5 (i) ;
     denu_x1(i)   =  numu1(i)*X(i*choices-4,2)  +  numu2(i)*X(i*choices-3,2) + numu3(i)*X(i*choices-2,2)+ ... 
                     numu4(i)*X(i*choices-1,2) +   numu5(i)*X(i*choices  ,2 ) ;
     prob_u1(i)     =  numu1 (i) /  denu(i)  ;     prob_u2(i)     =  numu2 (i) /  denu(i)  ;
     prob_u3(i)     =  numu3 (i) /  denu(i)  ;     prob_u4(i)     =  numu4 (i) /  denu(i)  ;
     prob_u5(i)     =  numu5 (i) /  denu(i)  ;  

     
     numf1 (i) =  exp(X(i*choices-4,3)*bbeta(7)    ) ;
     numf2 (i) =  exp(X(i*choices-3,3)*bbeta(7)  + 1*bbeta(1)  ) ;
     numf3 (i) =  exp(X(i*choices-2,3)*bbeta(7)  + 1*bbeta(2)  ) ;
     numf4 (i) =  exp(X(i*choices-1,3)*bbeta(7)  + 1*bbeta(3)  ) ;
     numf5 (i) =  exp(X(i*choices  ,3)*bbeta(7)  + 1*bbeta(4)  ) ;
     denf(i)   =  numf1 (i)  +  numf2 (i) + numf3 (i)+ numf4 (i)+ numf5 (i) ;
     denf_x1(i)   =  numf1(i)*X(i*choices-4,3)  +  numf2(i)*X(i*choices-3,3) + numf3(i)*X(i*choices-2,3)+ ... 
                     numf4(i)*X(i*choices-1,3) +   numf5(i)*X(i*choices  ,3 ) ;
     prob_f1(i)     =  numf1 (i) /  denf(i)  ;     prob_f2(i)     =  numf2 (i) /  denf(i)  ;
     prob_f3(i)     =  numf3 (i) /  denf(i)  ;     prob_f4(i)     =  numf4 (i) /  denf(i)  ;
     prob_f5(i)     =  numf5 (i) /  denf(i)  ; 

end
 
for i=1:(size(DATA ,1))/choices     
        
          if y_short(i)==1
          
           % score(i,1) =   - ze(i)*( prob_e2(i) )* X(i*choices,4 )- zu(i)*( prob_u2(i) )*1 - zf(i)*( prob_f2(i) )*1 ;
           score(i,1) =   - ze(i)*( prob_e2(i) )*X(i*choices,4)- zu(i)*( prob_u2(i) )*X(i*choices,4) - zf(i)*( prob_f2(i) )*X(i*choices,4) ;
           score(i,2) =   - ze(i)*( prob_e3(i) )*X(i*choices,5)- zu(i)*( prob_u3(i) )*X(i*choices,5) - zf(i)*( prob_f3(i) )*X(i*choices,5) ; 
           score(i,3) =   - ze(i)*( prob_e4(i) )*X(i*choices,6)- zu(i)*( prob_u4(i) )*X(i*choices,6) - zf(i)*( prob_f4(i) )*X(i*choices,6) ;
           score(i,4) =   - ze(i)*( prob_e5(i) )*X(i*choices,7)- zu(i)*( prob_u5(i) )*X(i*choices,7) - zf(i)*( prob_f5(i) )*X(i*choices,7) ;
            
           score(i,5) =     ze(i,1)*(X(i*choices-4,1)- (dene_x1(i) / dene(i)));
           score(i,6) =     zu(i,1)*(X(i*choices-4,2)- (denu_x1(i) / denu(i)));  
           score(i,7) =     zf(i,1)*(X(i*choices-4,3)- (denf_x1(i) / denf(i)));
           
          elseif y_short(i)==2
  
           score(i,1) =     ze(i)*(1-prob_e2(i))*X(i*choices,4) + zu(i)*(1-prob_u2(i))*X(i*choices,4) + zf(i)*(1- prob_f2(i))*X(i*choices,4) ;
           
           score(i,2) =   - ze(i)*( prob_e3(i) )*X(i*choices,5)- zu(i)*( prob_u3(i) )*X(i*choices,5) - zf(i)*( prob_f3(i) )*X(i*choices,5) ; 
           score(i,3) =   - ze(i)*( prob_e4(i) )*X(i*choices,6)- zu(i)*( prob_u4(i) )*X(i*choices,6) - zf(i)*( prob_f4(i) )*X(i*choices,6) ;
           score(i,4) =   - ze(i)*( prob_e5(i) )*X(i*choices,7)- zu(i)*( prob_u5(i) )*X(i*choices,7) - zf(i)*( prob_f5(i) )*X(i*choices,7) ;
             
           score(i,5) =     ze(i,1)*(X(i*choices-3,1)- (dene_x1(i) / dene(i)));
           score(i,6) =     zu(i,1)*(X(i*choices-3,2)- (denu_x1(i) / denu(i)));  
           score(i,7) =     zf(i,1)*(X(i*choices-3,3)- (denf_x1(i) / denf(i)));
           
           elseif y_short(i)==3
  
           score(i,2) =     ze(i)*(1-prob_e3(i))*X(i*choices,5) + zu(i)*(1-prob_u3(i))*X(i*choices,5) + zf(i)*(1- prob_f3(i))*X(i*choices,5) ;
           
           score(i,1) =   - ze(i)*( prob_e2(i) )*X(i*choices,4)- zu(i)*( prob_u2(i) )*X(i*choices,4) - zf(i)*( prob_f2(i) )*X(i*choices,4) ; 
           score(i,3) =   - ze(i)*( prob_e4(i) )*X(i*choices,6)- zu(i)*( prob_u4(i) )*X(i*choices,6) - zf(i)*( prob_f4(i) )*X(i*choices,6) ;
           score(i,4) =   - ze(i)*( prob_e5(i) )*X(i*choices,7)- zu(i)*( prob_u5(i) )*X(i*choices,7) - zf(i)*( prob_f5(i) )*X(i*choices,7) ;
             
           score(i,5) =     ze(i,1)*(X(i*choices-2,1)- (dene_x1(i) / dene(i)));
           score(i,6) =     zu(i,1)*(X(i*choices-2,2)- (denu_x1(i) / denu(i)));  
           score(i,7) =     zf(i,1)*(X(i*choices-2,3)- (denf_x1(i) / denf(i)));
           
           elseif y_short(i)==4
  
           score(i,3) =     ze(i)*(1-prob_e4(i))*X(i*choices,6) + zu(i)*(1-prob_u4(i))*X(i*choices,6) + zf(i)*(1- prob_f4(i))*X(i*choices,6) ;
           
           score(i,1) =   - ze(i)*( prob_e2(i) )*X(i*choices,4)- zu(i)*( prob_u2(i) )*X(i*choices,4) - zf(i)*( prob_f2(i) )*X(i*choices,4) ; 
           score(i,2) =   - ze(i)*( prob_e3(i) )*X(i*choices,5)- zu(i)*( prob_u3(i) )*X(i*choices,5) - zf(i)*( prob_f3(i) )*X(i*choices,5) ;
           score(i,4) =   - ze(i)*( prob_e5(i) )*X(i*choices,7)- zu(i)*( prob_u5(i) )*X(i*choices,7) - zf(i)*( prob_f5(i) )*X(i*choices,7) ;
             
           score(i,5) =     ze(i,1)*(X(i*choices-1,1)- (dene_x1(i) / dene(i)));
           score(i,6) =     zu(i,1)*(X(i*choices-1,2)- (denu_x1(i) / denu(i)));  
           score(i,7) =     zf(i,1)*(X(i*choices-1,3)- (denf_x1(i) / denf(i)));

           elseif y_short(i)==5
  
           score(i,4) =     ze(i)*(1-prob_e5(i))*X(i*choices,7) + zu(i)*(1-prob_u5(i))*X(i*choices,7) + zf(i)*(1- prob_f5(i))*X(i*choices,7) ;
           
           score(i,1) =   - ze(i)*( prob_e2(i) )*X(i*choices,4)- zu(i)*( prob_u2(i) )*X(i*choices,4) - zf(i)*( prob_f2(i) )*X(i*choices,4) ; 
           score(i,2) =   - ze(i)*( prob_e3(i) )*X(i*choices,5)- zu(i)*( prob_u3(i) )*X(i*choices,5) - zf(i)*( prob_f3(i) )*X(i*choices,5) ;
           score(i,3) =   - ze(i)*( prob_e4(i) )*X(i*choices,6)- zu(i)*( prob_u4(i) )*X(i*choices,6) - zf(i)*( prob_f4(i) )*X(i*choices,6) ;
             
           score(i,5) =     ze(i,1)*(X(i*choices ,1)- (dene_x1(i) / dene(i)));
           score(i,6) =     zu(i,1)*(X(i*choices ,2)- (denu_x1(i) / denu(i)));  
           score(i,7) =     zf(i,1)*(X(i*choices ,3)- (denf_x1(i) / denf(i)));
                     
          else
              disp 'not good'
              pause
          end
     
end
  
sscore = sum(score) ;
loglik = [sscore ];
  
  
  
  